

# CODE TO MAKE PLOTS FOR FIGURE 4
# 	This code will make each panel in Figure 4, which we then combined in Illustrator to make the final figure

rm(list=ls())
library(class)  #user will need to install class package if not already installed
wd <- ""  #insert into quotation marks the working directory to location where replication file was unzipped. E.g /Documents/BDLMS_Replication/ 	Make sure to include the final backslash.
"%&%"<-function(x,y)paste(x,y,sep="")  #define a function for easy string pasting


# Again we are going to fill some lists that we will use to generate values in the final columns of Table 2
numstud=7
tail=vector("list",numstud)  #going to collect tails in this list. tails of 95% CI, not full tails.
	# for study i, tail[[i]][1:3] gives negative tails for reg uncer, clim uncer, total uncer
	#		and tail[[i]][4:6] gives positive tail for same.  
studnms=vector("list",numstud)  #collect study names

############################################################
# 	Dell Jones Olken paper (DJO) - AEJ Macro 2012
############################################################

# Unlike for the other figures, we actually calculate these in a separate script since they are a little more complicated. Script is "djo_projections.R"

inlist=1

setwd(wd%&%"data/DJO/")  #navigate to DJO data directory


load("DJO_projections.Rdata")  #the impact projections we calculated
boot=read.csv("boot_djo.csv")  #boostrapped reg coefficients from DJO specification in Table 2, Col 3
boot=boot[,2] #just keep temp coefficients

#collect scenario and model names
nms=strsplit(dimnames(impact)[[2]],"\\_")
scen=mod=per=c()
for (i in 1:length(nms)) {
	scen=c(scen,nms[[i]][1]) #scenario
	mod=c(mod,nms[[i]][2]) #climate model
	per=c(per,nms[[i]][3]) #mid- or end-of-century
}

# location of median bootstrap temperature coefficient in impact array
medloc=as.numeric(knn1(boot,median(boot),1:1000))

pdf(file=wd%&%"output/Fig4subfigs/DJO_uncertainty_plot.pdf",width=5,height=5)
par(mar=c(4,3,2,1))
plot(1,type="n",xlim=c(-90,0),ylim=c(0.5,8),yaxt="n",xlab="change in poor country GDP/cap (%)",ylab="",cex.axis=0.8, main="DJO, 2012",xaxt="n")
sq=seq(-90,0,10)
axis(1,at=sq,sq,cex.axis=0.8)
abline(v=0,lwd=2,col="black",lty=3)
abline(v=sq[sq!=0],lwd=1,col="grey",lty=3)
rect(-100,3.65,10,4.35,col="black")
text(-45,4,"2080-2100",cex=1.2,col="white")
rect(-100,7.65,10,8.35,col="black")
text(-45,7.95,"2040-2060",cex=1.2,col="white")
box()
mat=matrix(c(5:7,1:3),byrow=T,nrow=2)

# Loop over time periods (mid- or end-century)
pers=c("mid","end")
for (p in 1:2) {

uncer=c()  #list that we will fill
scenper=which(scen=="a1b"&per==pers[p])  #pulls the models in the scenario and time period we care about
med=impact[medloc,scenper]

a1b_med=which(med==sort(med)[length(med)/2])  #which a1b model is the median when evaluated at the regression coefficients
uncer[[1]]=impact[,scenper[a1b_med]]  #bootstraps evaluated at the median a1b model
uncer[[2]]=med  #climate model uncertainty
uncer[[3]]=c(impact[,scenper])  #bootstraps evaluated at all the a1b models

colz=c("grey99","grey75","grey50")

rg=c()
for (n in 1:3) {
	toplot=c(uncer[[n]])
	low=quantile(toplot,probs=c(0.025)) #to make whiskers right for 95% conf interval, we reset the tails to the values of the 2.5% and 97.5%.
	hi=quantile(toplot,probs=c(0.975))
	toplot=replace(toplot,toplot<low,low)
	toplot=replace(toplot,toplot>hi,hi)
	boxplot(toplot,horizontal=T,at=mat[p,n],add=T,range=0,outline=FALSE,axes=F,col=colz[n])
	rg=c(rg,max(toplot)-min(toplot))
	if (p==1) {tail[[inlist]][c(n,n+3)]=range(toplot)}
	}
	text(-90,mat[p,1],paste("[",round(rg[3]/rg[1],2)*100-100,"%]",sep=""),pos=4,cex=1)

}
dev.off()
studnms[[inlist]]="DJO"

####################################################################
# 	Mendelsohn Nordhaus Shaw (MNS), AER 1994
#	and Schlenker Hanemann Fischer (SHF), AER 2005
####################################################################

setwd(wd%&%"data/MNS_SHF/")  #navigate to MNS/SHF data directory

# Bring in the data: boostrapped regression coefficients created in "MNS_uncertainty.do", and country level climate projections that are weighted by crop area (created in get_climchg_MNS.R)

inlist=2
mns=read.csv("boot_mns.csv")  #bootstrapped regression coefficients, column 2 in Table 3 of MNS
shf=read.csv("boot_shf.csv")  # same specification, but restricted to non-irrigated non-urban counties as in SHF
data=read.csv("mns_sample.csv")
load("MNS_SHF_climchg_temp.Rdata")
tm=tm*9/5  # MNS and SHF temperature variables are in farenheight, but climate change projections are in C
load("MNS_SHF_climchg_prec.Rdata")
baseprec=data.frame(data$janprec_raw,data$aprprec_raw,data$julprec_raw,data$octprec_raw)

# GENERATE PLOT FOR MNS RESULTS
pdf(file=wd%&%"output/Fig4subfigs/MNS_uncertainty_plot.pdf",width=5,height=5)
par(mar=c(4,3,2,1))
plot(1,type="n",xlim=c(-800,1100),ylim=c(0.5,8),yaxt="n",xaxt="n",xlab="change in land values ($bil)",ylab="",cex.axis=0.8, main="MNS, 1994")
sq=seq(-800,1000,200)
axis(1,at=sq,sq,cex.axis=0.7)
abline(v=0,lwd=2,col="black",lty=3)
abline(v=sq[sq!=0],lwd=1,col="grey",lty=3)
rect(-1000,3.65,1200,4.35,col="black")
text(median(sq),4,"2080-2100",cex=1.2,col="white")
rect(-1000,7.65,1200,8.35,col="black")
text(median(sq),7.95,"2040-2060",cex=1.2,col="white")
box()
mat=matrix(c(5:7,1:3),byrow=T,nrow=2)
notna=(is.na(data$cropsales)==F)

#loop over time periods
for (p in 1:2) {

# construct values at which to evaluate the boootstrapped coefficients
m1=as.matrix(mns[,2:17])*sum(data$farmland[notna])/1e9  #boostrapped estimates for the 16 independent variables, converted to billions of dollars (regression coeffs are in dollars/acre, so multiplying by total acres and dividing by a billion)
m2=c()  #the values at which to evaluate the boootstrapped estimates
for (i in 1:4) {
	m2=rbind(m2,tm[,i,p,2],tm[,i,p,2]^2)
	}
for (i in 1:4) {
	prchg=pr[,i,p,2]*weighted.mean(baseprec[notna,i],data$farmland[notna])  #projections are in %chg, converting to levels to evaluate coeffs
	m2=rbind(m2,prchg,prchg^2)
	}

mods=dimnames(tm)[[1]]  #modelnames
a1b=grep("a1b",mods)  # which of these models are a1b
uncer=c()
# get medians
medboot=apply(m1,2,median)  #median bootstrapped coefficients, about the same as reg coefficients
med=medboot%*%m2[,a1b]  #climate uncertainty, evaluated at median bootstrapped coefficients
a1b_med=which(med==sort(med)[length(med)/2])  #which a1b model is the median when evaluated at the regression coefficients
uncer[[1]]=m1%*%m2[,a1b_med]  #boostraps evaluated at the median a1b model
uncer[[2]]=med
uncer[[3]]=m1%*%m2[,a1b]  #bootstraps evaluated at all the a1b models

# now plot boxplots
bottom=sum(as.numeric(data$farmland1982[notna]))*weighted.mean(data$farmval1982[notna],data$cropsales[notna])/1e9  #gives current value of land in billions.  using this as lower bound so losses can exceed 100%
colz=c("grey99","grey75","grey50")
rg=c()
for (n in 1:3) {
	toplot=c(uncer[[n]])
	low=quantile(toplot,probs=c(0.025)) #to make whiskers right for 95% conf interval, we reset the tails to the values of the 2.5% and 97.5%.
	hi=quantile(toplot,probs=c(0.975))
	toplot=replace(toplot,toplot<low,low)
	toplot=replace(toplot,toplot>hi,hi)
	toplot=replace(toplot,toplot<(-bottom),-bottom)  #constrain losses to existing value of farmland
	boxplot(toplot,horizontal=T,at=mat[p,n],add=T,range=0,outline=FALSE,axes=F,col=colz[n])
	rg=c(rg,max(toplot)-min(toplot))
	if (p==1) {tail[[inlist]][c(n,n+3)]=range(toplot)}
	}
	text(-850,mat[p,1],paste("[",round(rg[3]/rg[1],2)*100-100,"%]",sep=""),pos=4,cex=1)

}
dev.off()
studnms[[inlist]]="MNS"

# GENERATE PLOT FOR SHF RESULTS
inlist=3
pdf(file=wd%&%"output/Fig4subfigs/SHF_uncertainty_plot.pdf",width=5,height=5)
par(mar=c(4,3,2,1))
plot(1,type="n",xlim=c(-600,130),ylim=c(0.5,8),yaxt="n",xaxt="n",xlab="change in land values ($bil)",ylab="",cex.axis=0.8, main="SHF, 2005")
sq=seq(-600,100,100)
axis(1,at=sq,sq,cex.axis=0.7)
abline(v=0,lwd=2,col="black",lty=3)
abline(v=sq[sq!=0],lwd=1,col="grey",lty=3)
rect(-700,3.65,400,4.35,col="black")
text(-200,4,"2080-2100",cex=1.2,col="white")
rect(-700,7.65,400,8.35,col="black")
text(-200,7.95,"2040-2060",cex=1.2,col="white")
box()
mat=matrix(c(5:7,1:3),byrow=T,nrow=2)

smpl=data$dryrural==1

#loop over time periods
for (p in 1:2) {

# construct values at which to evaluate the boootstrapped coefficients
m1=as.matrix(shf[,2:17])*sum(data$farmland[smpl])/1e9  #boostrapped estimates for the 16 independent variables, converted to billions of dollars (regression coeffs are in dollars/acre, so multiplying by total acres and dividing by a billion)
m2=c()  #the values at which to evaluate the boootstrapped estimates
for (i in 1:4) {
	m2=rbind(m2,tm[,i,p,2],tm[,i,p,2]^2)
	}
for (i in 1:4) {
	prchg=pr[,i,p,2]*weighted.mean(baseprec[smpl,i],data$farmland[smpl])  #projections are in %chg, converting to levels to evaluate coeffs
	m2=rbind(m2,prchg,prchg^2)
	}

mods=dimnames(tm)[[1]]  #modelnames
a1b=grep("a1b",mods)  # which of these models are a1b
uncer=c()
# get medians
medboot=apply(m1,2,median)  #median bootstrapped coefficients, about the same as reg coefficients
med=medboot%*%m2[,a1b]  #climate uncertainty, evaluated at median bootstrapped coefficients
a1b_med=which(med==sort(med)[length(med)/2])  #which a1b model is the median when evaluated at the regression coefficients
uncer[[1]]=m1%*%m2[,a1b_med]  #boostraps evaluated at the median a1b model
uncer[[2]]=med
uncer[[3]]=m1%*%m2[,a1b]  #bootstraps evaluated at all the a1b models

bottom=sum(as.numeric(data$farmland1982[smpl]))*weighted.mean(data$farmval1982[smpl],data$totalcropland1982[smpl])/1e9  #gives current value of land in billions.  using this as lower bound so losses can exceed 100%
colz=c("grey99","grey75","grey50")
rg=c()
for (n in 1:3) {
	toplot=c(uncer[[n]])
	low=quantile(toplot,probs=c(0.025)) #to make whiskers right for 95% conf interval, we reset the tails to the values of the 2.5% and 97.5%.
	hi=quantile(toplot,probs=c(0.975))
	toplot=replace(toplot,toplot<low,low)
#	toplot=replace(toplot,toplot<(-100),(-100))  #use this line if using % changes
	toplot=replace(toplot,toplot<(-bottom),-bottom)  #use this line if using changes in billions
	toplot=replace(toplot,toplot>hi,hi)
	boxplot(toplot,horizontal=T,at=mat[p,n],add=T,range=0,outline=FALSE,axes=F,col=colz[n])
	rg=c(rg,max(toplot)-min(toplot))
	if (p==1) {tail[[inlist]][c(n,n+3)]=range(toplot)}
	}
	text(-630,mat[p,1],paste("[",round(rg[3]/rg[1],2)*100-100,"%]",sep=""),pos=4,cex=1)

}
dev.off()
studnms[[inlist]]="SHF"



####################################################################
# 	DESCHENES AND GREENSTONE 2007 AER
####################################################################

inlist=4
setwd(wd%&%"data/DG/")  #navigate to DG data directory

load("DG_climchg.Rdata")
boot=read.csv("boot_dg.csv")
data=read.csv("dg_data.csv")

m1=as.matrix(boot[,2:9])
#coefficients are in $/acre, and want to convert to overall damages in $bil, so need to multiply by total acres.
# multiplying dryland coefficients by dryland acres, and irrigated coefficients by irrigated acres, and dividing by 1 billion
tot=aggregate(data,by=list(data$fips),mean,na.rm=T)
dryacre=sum(tot$fland[tot$dry==1])
irracre=sum(tot$fland[tot$dry==0])
m1[,1:4]=m1[,1:4]*dryacre/1e9
m1[,5:8]=m1[,5:8]*irracre/1e9

mods=dimnames(climchg)[[1]]  #modelnames
a1b=grep("a1b",mods)  # which of these models are a1b
medboot=apply(m1,2,median)  #median bootstrapped coefficients, about the same as reg coefficients

#make plot
pdf(file=wd%&%"output/Fig4subfigs/DG_uncertainty_plot.pdf",width=5,height=5)
par(mar=c(4,3,2,1))
plot(1,type="n",xlim=c(-5,4),ylim=c(0.5,8),yaxt="n",xaxt="n",xlab="change in profits ($bil)",ylab="",cex.axis=0.8, main="DG, 2007")
sq=seq(-10,6,2)
axis(1,at=sq,sq,cex.axis=0.85)
abline(v=0,lwd=2,col="black",lty=3)
abline(v=sq[sq!=0],lwd=1,col="grey",lty=3)
rect(-12,3.65,10,4.35,col="black")
text(-0.7,4,"2080-2100",cex=1.2,col="white")
rect(-12,7.65,10,8.35,col="black")
text(-0.7,7.95,"2040-2060",cex=1.2,col="white")
box()
mat=matrix(c(5:7,1:3),byrow=T,nrow=2)

#loop over periods
for (p in 1:2) {

uncer=c()
m2=t(climchg[,,p])
med=medboot%*%m2[,a1b]  #climate uncertainty, evaluated at median bootstrapped coefficients
a1b_med=which(med==sort(med)[length(med)/2])  #which a1b model is the median when evaluated at the regression coefficients
uncer[[1]]=m1%*%m2[,a1b_med]  #boostraps evaluated at the median a1b model
uncer[[2]]=med
uncer[[3]]=m1%*%m2[,a1b]  #bootstraps evaluated at all the a1b models

colz=c("grey99","grey75","grey50")
rg=c()
for (n in 1:3) {
	toplot=c(uncer[[n]])
	low=quantile(toplot,probs=c(0.025)) #to make whiskers right for 95% conf interval, we reset the tails to the values of the 2.5% and 97.5%.
	hi=quantile(toplot,probs=c(0.975))
	toplot=replace(toplot,toplot<low,low)
#	toplot=replace(toplot,toplot<(-bottom),-bottom)  #use this line if using changes in billions
	toplot=replace(toplot,toplot>hi,hi)
	boxplot(toplot,horizontal=T,at=mat[p,n],add=T,range=0,outline=FALSE,axes=F,col=colz[n])
	rg=c(rg,max(toplot)-min(toplot))
	if (p==1) {tail[[inlist]][c(n,n+3)]=range(toplot)}
	}
	text(-5.2,mat[p,1],paste("[",round(rg[3]/rg[1],2)*100-100,"%]",sep=""),pos=4,cex=1)
}
dev.off()
studnms[[inlist]]="DG"


####################################################################
# 	FISHER HANEMANN ROBERTS SCHLENKER 2010 AER COMMENT
####################################################################

inlist=5

setwd(wd%&%"data/FHRS/")  #navigate to DJO data directory

load("FHRS_climchg.Rdata")
boot=read.csv("boot_fhrs.csv")
data=read.csv("fhrs_data.csv")

m1=as.matrix(boot[,2:9])
#coefficients are in bushels/acre, and want to convert to overall % change in yields.  going to do this by converting into % effect on overall production across all counties, which will then allow us to easily sum across dryland and irrigated counties
#average across years
tot=aggregate(data,by=list(data$fips),mean,na.rm=T)
totprod=sum(tot$corn_planted*tot$yield,na.rm=T)

# multiply area-weighted yield effect (what the regression spits out) by corn area to get effect of change on overall production, then lower down we divide by total production
m1[,1:4]=m1[,1:4]*sum(tot$corn_planted[tot$dry==1],na.rm=T)/totprod*100
m1[,5:8]=m1[,5:8]*sum(tot$corn_planted[tot$dry==0],na.rm=T)/totprod*100

mods=dimnames(climchg)[[1]]  #modelnames
a1b=grep("a1b",mods)  # which of these models are a1b
medboot=apply(m1,2,median)  #median bootstrapped coefficients, about the same as reg coefficients

#make plot
pdf(file=wd%&%"output/Fig4subfigs/FHRS_uncertainty_plot.pdf",width=5,height=5)
par(mar=c(4,3,2,1))
plot(1,type="n",xlim=c(-80,0),ylim=c(0.5,8),yaxt="n",xaxt="n",xlab="change in corn yields (%)",ylab="",cex.axis=0.8, main="FHRS, 2010")
sq=seq(-100,0,20)
axis(1,at=sq,sq,cex.axis=0.85)
abline(v=0,lwd=2,col="black",lty=3)
abline(v=sq[sq!=0],lwd=1,col="grey",lty=3)
rect(-120,3.65,10,4.35,col="black")
text(-40,4,"2080-2100",cex=1.2,col="white")
rect(-120,7.65,10,8.35,col="black")
text(-40,7.95,"2040-2060",cex=1.2,col="white")
box()
mat=matrix(c(5:7,1:3),byrow=T,nrow=2)

#loop over periods
for (p in 1:2) {

uncer=c()
m2=t(climchg[,,p])
med=medboot%*%m2[,a1b]  #climate uncertainty, evaluated at median bootstrapped coefficients
a1b_med=which(med==sort(med)[length(med)/2])  #which a1b model is the median when evaluated at the regression coefficients
uncer[[1]]=m1%*%m2[,a1b_med]  #boostraps evaluated at the median a1b model
uncer[[2]]=med
uncer[[3]]=m1%*%m2[,a1b]  #bootstraps evaluated at all the a1b models

colz=c("grey99","grey75","grey50")
rg=c()
for (n in 1:3) {
	toplot=c(uncer[[n]])
	low=quantile(toplot,probs=c(0.025)) #to make whiskers right for 95% conf interval, we reset the tails to the values of the 2.5% and 97.5%.
	hi=quantile(toplot,probs=c(0.975))
	toplot=replace(toplot,toplot<low,low)
	toplot=replace(toplot,toplot<(-100),-100)  
	toplot=replace(toplot,toplot>hi,hi)
	boxplot(toplot,horizontal=T,at=mat[p,n],add=T,range=0,outline=FALSE,axes=F,col=colz[n])
	rg=c(rg,max(toplot)-min(toplot))
	if (p==1) {tail[[inlist]][c(n,n+3)]=range(toplot)}
	}
	text(-80,mat[p,1],paste("[",round(rg[3]/rg[1],2)*100-100,"%]",sep=""),pos=4,cex=1)
}
dev.off()
studnms[[inlist]]="FHRS"



############################################################
# 	Burke et al (BMSDL) - PNAS 2009
############################################################

inlist=6
setwd(wd%&%"data/BMSDL/")  #navigate to BMSDL data directory

boot=read.csv("boot_BMSDL.csv")  #boostrapped reg coefficients from pnas specification 1 in Table 1 of PNAS paper
temp=read.csv("BMSDL_climchg.csv")  #temperature projections

# convert coefficients to %change
boot=boot[,2]+boot[,3]  #evaluating them at same temp change, so can just add together first to simplify
boot=boot/0.11*100  #0.11 is the baseline incidence of conflict in the sample

# Generate temperature projections at which to evaluate regression coefficients. Regressions are unweighted, so temp projections are just a simple average over projections for countires in the sample. 
tchg=apply(temp[,4:dim(temp)[2]],2,mean)

#collect scenario and model names
nms=strsplit(names(tchg),"\\_")
scen=mod=per=c()
for (i in 1:length(nms)) {
	scen=c(scen,nms[[i]][1]) #scenario
	mod=c(mod,nms[[i]][2]) #climate model
	per=c(per,nms[[i]][3]) #mid- or end-of-century
}

# median bootstrap coefficients
m1=matrix(boot,nrow=length(boot))
medboot = apply(m1,2,median)

pdf(file=wd%&%"output/Fig4subfigs/BMSDL_uncertainty_plot.pdf",width=5,height=5)
par(mar=c(4,3,2,1))
plot(1,type="n",xlim=c(-70,380),ylim=c(0.5,8),yaxt="n",xlab="increase in conflict incidence (%)",ylab="",cex.axis=0.8, main="BMSDL, 2009",xaxt="n")
sq=seq(-50,350,50)
axis(1,at=sq,sq,cex.axis=0.8)
abline(v=0,lwd=2,col="black",lty=3)
abline(v=sq[sq!=0],lwd=1,col="grey",lty=3)
rect(-100,3.65,450,4.35,col="black")
text(150,4,"2080-2100",cex=1.2,col="white")
rect(-100,7.65,450,8.35,col="black")
text(150,7.95,"2040-2060",cex=1.2,col="white")
box()
mat=matrix(c(5:7,1:3),byrow=T,nrow=2)

# Loop over time periods (mid- or end-century)
pers=c("mid","end")
for (p in 1:2) {

uncer=c()  #list that we will fill
scenper=which(scen=="a1b"&per==pers[p])  #pulls the models in the scenario and time period we care about

# climate model uncertainty for a1b and period in question
m2=matrix(tchg,ncol=length(tchg))  #temperature changes we're evaluating
med=medboot*m2[,scenper]  #climate uncertainty across the models and time period in question

a1b_med=which(med==sort(med)[length(med)/2])  #which a1b model is the median when evaluated at the regression coefficients
uncer[[1]]=m1%*%m2[,scenper[a1b_med]]  #bootstraps evaluated at the median a1b model
uncer[[2]]=med
uncer[[3]]=m1%*%m2[,scenper]  #bootstraps evaluated at all the a1b models

colz=c("grey99","grey75","grey50")

rg=c()
for (n in 1:3) {
	toplot=c(uncer[[n]])
	low=quantile(toplot,probs=c(0.025)) #to make whiskers right for 95% conf interval, we reset the tails to the values of the 2.5% and 97.5%.
	hi=quantile(toplot,probs=c(0.975))
	toplot=replace(toplot,toplot<low,low)
	toplot=replace(toplot,toplot>hi,hi)
	boxplot(toplot,horizontal=T,at=mat[p,n],add=T,range=0,outline=FALSE,axes=F,col=colz[n])
	rg=c(rg,max(toplot)-min(toplot))
	if (p==1) {tail[[inlist]][c(n,n+3)]=range(toplot)}
	}
	text(-90,mat[p,1],paste("[",round(rg[3]/rg[1],2)*100-100,"%]",sep=""),pos=4,cex=1)

}
dev.off()
studnms[[inlist]]="BMSDL"


############################################################
# 	Schlenker and Lobell ERL 2010 
############################################################

inlist=7
setwd(wd%&%"data/SL/")  #navigate to SL data directory

clim=read.csv('SL_climchg.csv')
boot=read.csv('boot_sl.csv')
data=read.csv('africa_yield_clim.csv')

prec_base=aggregate(data$maize_prec,by=list(data$fidcode),mean,na.rm=T)  #country-average precip
prec_base=mean(prec_base[,2])

tchg=apply(clim[,2:103],2,mean)  #regressions are unweighted so evaluating at a simple average of temp and precip change across countries
pchg=apply(clim[,104:dim(clim)[2]],2,mean)*prec_base  #precip changes are in %, so evaluating at continental mean annual precip to get changes in levels. 

#collect scenario and model names
nms=strsplit(names(tchg),"\\_")
scen=mod=per=c()
for (i in 1:length(nms)) {
	scen=c(scen,nms[[i]][2]) #scenario
	mod=c(mod,nms[[i]][3]) #climate model
	per=c(per,nms[[i]][4]) #mid- or end-of-century
}

# median bootstrap coefficients
m1=data.matrix(boot[,2:3])
medboot = apply(m1,2,median)

#climate change projections
m2=data.matrix(rbind(tchg,pchg))

pdf(file=wd%&%"output/Fig4subfigs/SL_uncertainty_plot.pdf",width=5,height=5)
par(mar=c(4,3,2,1))
plot(1,type="n",xlim=c(-60,0),ylim=c(0.5,8),yaxt="n",xlab="change in corn yield (%)",ylab="",cex.axis=0.8, main="SL, 2010",xaxt="n")
sq=seq(-60,0,10)
axis(1,at=sq,sq,cex.axis=0.8)
abline(v=0,lwd=2,col="black",lty=3)
abline(v=sq[sq!=0],lwd=1,col="grey",lty=3)
rect(-100,3.65,20,4.35,col="black")
text(-30,4,"2080-2100",cex=1.2,col="white")
rect(-100,7.65,20,8.35,col="black")
text(-30,7.95,"2040-2060",cex=1.2,col="white")
box()
mat=matrix(c(5:7,1:3),byrow=T,nrow=2)

# Loop over time periods (mid- or end-century)
pers=c("mid","end")
for (p in 1:2) {

uncer=c()  #list that we will fill
scenper=which(scen=="a1b"&per==pers[p])  #pulls the models in the scenario and time period we care about

# climate model uncertainty for a1b and period in question
med=medboot%*%m2[,scenper]  #climate uncertainty across the models and time period in question

a1b_med=which(med==sort(med)[length(med)/2])  #which a1b model is the median when evaluated at the regression coefficients
uncer[[1]]=m1%*%m2[,scenper[a1b_med]]  #bootstraps evaluated at the median a1b model
uncer[[2]]=med
uncer[[3]]=m1%*%m2[,scenper]  #bootstraps evaluated at all the a1b models

colz=c("grey99","grey75","grey50")

rg=c()
for (n in 1:3) {
	toplot=c(uncer[[n]])*100
	low=quantile(toplot,probs=c(0.025)) #to make whiskers right for 95% conf interval, we reset the tails to the values of the 2.5% and 97.5%.
	hi=quantile(toplot,probs=c(0.975))
	toplot=replace(toplot,toplot<low,low)
	toplot=replace(toplot,toplot>hi,hi)
	boxplot(toplot,horizontal=T,at=mat[p,n],add=T,range=0,outline=FALSE,axes=F,col=colz[n])
	rg=c(rg,max(toplot)-min(toplot))
	if (p==1) {tail[[inlist]][c(n,n+3)]=range(toplot)}
	}
	text(-60,mat[p,1],paste("[",round(rg[3]/rg[1],2)*100-100,"%]",sep=""),pos=4,cex=1)

}
dev.off()
studnms[[inlist]]="SL"


############################
#  Analysis of tails for 2050
#############################

print("Factor by which lower tail of total uncertainty exceeds lower tail of regression uncertainty, A1b, 2050")
for (i in 1:numstud) {
	rg=tail[[i]]
	print(paste(studnms[[i]],": ",round(rg[3]/rg[1],1),sep=""))
}
